Thermodynamics of a bad metal-Mott insulator transition in the presence of frustration 
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We study a range of thermodynamic properties (charge susceptibility, specific heat, entropy and spin sus- 
ceptibility) of the Hubbard model on the anisotropic triangular lattice at half filling by means of the numerical 
finite-temperature Lanczos method. We observe clear signatures of a metal-Mott insulator transition in the 
charge susceptibility, and show that the metallic phase is characterized by a small charge susceptibility, large 
entropy, low coherence temperature, large renormalized quasiparticle mass, and large spin susceptibility. We 
show that the local magnetic moment in the metallic phase is large and comparable to the local moment in the 
insulating phase. These are all characteristics of bad metallic behavior. In addition, we show how frustration 
increases the density of low-lying spin excitations in the Mott insulating phase and decreases longer range spin 
correlations. 

PACS numbers: 71.72.+a, 71.30.+h, 74.25.-q, 74.70.Kn, 75.20.-g 



Remarkable observations of a possible spin liquid phase [1] 
and a new universality class of the metal-insulator transition 
[2] in organic charge transfer salts, which in addition show 
unconventional superconductivity [3], have increased interest 
in these materials. It has been argued that a proper micro- 
scopic description of these material can be given with a Hub- 
bard model on the anisotropic triangular lattice at half filling 
[4]. Parameters of the model for the description of organic 
charge transfer salts fall into the regime of strong correlations 
and significant frustration of antiferromagnetic spin interac- 
tions. This is the most challenging parameter regime, where 
analytical approaches become unreliable, and one needs to re- 
sort to numerical techniques. 

In this Letter we study a range of thermodynamic proper- 
ties (charge susceptibility, specific heat, entropy and spin sus- 
ceptibility) of the Hubbard model on the anisotropic triangu- 
lar lattice at half-filling. The model exhibits a Mott metal- 
insulator transition (MIT), which can be driven either by in- 
teraction strength or by frustration. We argue that the metal- 
lic phase has a strongly reduced coherence temperature T CO h, 
below which a Fermi liquid metal with coherent quasiparti- 
cle excitations may exist. Above T co h the model is in a bad 
metallic regime with large local magnetic moments. We show 
how frustration increases the low temperature specific heat, 
entropy and spin susceptibility in the insulating phase. Al- 
though the charge susceptibility shows definitive signatures of 
the metal-insulator transition, the specific heat and spin sus- 
ceptibility do not. Indeed, above T co h there appears to be little 
difference between the bad metal and the Mott insulator. This 
is similar to the dynamical mean-field theory (DMFT) picture 
of the transition [5, 6]. 

Model. The Hubbard model on the anisotropic triangular 
lattice has the Hamiltonian 

The hopping parameters tij = t for nearest neighbors in two 
directions of the triangular lattice, while tij = t' for nearest 
neighbors in the third direction. Cj jCT (cj ) is a fermionic anni- 
hilation (creation) operator for an electron on site i with spin 



a (either j" or 4-). = c] Ci t(T , U is the on-site Coulomb 
repulsion, and fj, is the chemical potential. Most of our results 
are presented in units of t, and we use h = ks = 1. We only 
consider the case of half-filling since this is relevant to several 
important families of organic charge transfer salts [4]. 

Numerical method. To calculate thermodynamic proper- 
ties for the model we use the finite-temperature Lanczos 
method (FTLM) [7-9]. Within FTLM, the Hamiltonian is 
effectively diagonalized on a small cluster. We use 16 site 
clusters with twisted boundary conditions. With FTLM one is 
essentially capable of obtaining results in the thermodynamic 
limit for temperatures T above some limiting temperature Tf s , 
below which finite size effects become important. Due to ge- 
ometrical frustration on the triangular lattice, the spin corre- 
lation length is substantially reduced [10], making finite size 
effects and Tf S smaller. Additionally, we apply averaging over 
twisted boundary conditions [11, 12], which further reduces 
finite size effects, and allows exact results to be obtained in 
the thermodynamic limit for U = [12]. Our results are ob- 
tained with a number of samples over different twisted bound- 
ary conditions Ng > 32, while the effect of averaging over 
random vectors used in the FTLM plays only a minor role. 
This is due to the large number of many-body states and finite 
T, and therefore usually one random vector suffices [8, 12]. 
More details on the calculation of thermodynamic properties 
with FTLM can be found in the Supplementary Material [13] 
and Refs. 8 and 12. Recently, it was shown that for the t-J 
model, the FTLM gives results in agreement with those ob- 
tained by a numerical linked-cluster algorithm [14] suggesting 
that FTLM on small lattices can give results comparable to the 
thermodynamic limit for strongly correlated metallic phases. 

Charge susceptibility. In Fig. 1 we show the temper- 
ature dependence of the charge susceptibility \ c = 
which is strongly suppressed with increasing U from its non- 
interacting electron value (calculated for infinite system). This 
is primarily due to broadening of the density of states over a 
larger energy range (W + U) or due to reduced quasiparticle 
weight [13]. In addition to this overall decrease of \c with in- 
creasing U, Xc becomes further suppressed at low T for higher 
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Figure 1. (color online) Signatures of the Mott transition in the 
temperature dependence of the charge susceptibility \ c . The figure 
shows Xc vs. temperature T for several interaction strengths U and 
t' — t. Xc decreases with increasing U and is almost independent 
of temperature in the metallic phase (U < It). In contrast, in the 
insulating phase (U > 8t) it is strongly suppressed at low T, with an 
activated behavior x^' = ae~ Ac//T (fits shown with thin red lines). 



gence for U — >• U c in the metallic regime (U < U c ). This is 
in contrast to what was observed for a filling-controlled MIT 
within DMFT [15] and a path-integral renormalization group 
approach on the square lattice with next-nearest-neighbor 
hopping [16, 17]. On the other hand, no sign of divergence 
was observed for filling controlled MIT in an exact diago- 
nalization study on the triangular lattice [18]. This suggests 
that upon changing from a filling- to bandwidth-controlled 
transition, dimensionality or frustration can affect the type of 
Mott MIT. The possibility of different characters of filling- 
and bandwidth-controlled MIT was pointed out in Ref. 16. 

Although our results at finite T do not allow precise deter- 
mination of the order of the MIT, the linear dependence of A c 
on U (Fig. 2), which persists down to U quite close to U c , 
is in agreement with a V- or T-shaped metal-insulator bound- 
ary in the fi-U plane and therefore also in agreement with the 
suggested [16, 19] first order transition. However, our results 
cannot rule out a second order phase transition, as proposed 
by Senfhil[20]. 
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Figure 2. (color online) Signatures of the Mott transition in the 
charge susceptibility Xc at T = O.lt and charge gap A c . Fig- 
ure shows the reduction in Xc(T = O.lt) with increasing U for 
several frustrations t'/t and the charge gap A c vs. U for t' — t. 
Xc(T = O.lt) is finite and slowly decreasing with increasing U at 
smaller values of U (< U c ), which corresponds to metallic behavior. 
At some larger value of U ~ U c , Xc(T = O.lt) becomes strongly 
suppressed exhibiting a metal-insulator transition (MIT) and is close 
to zero for larger U (> U c ), corresponding to Mott insulating behav- 
ior. In this regime a charge gap opens, which increases with increas- 
ing U. The figure also shows that the MIT appears at smaller U c for 
less frustrated systems (smaller t'/t). This is observed also with the 
move of A c curves to the left with decreasing t'/t (not shown). 



U > U c due to the MIT and opening of a charge gap A c (see 
Fig. 2). In the insulating phase, Xc shows an activated be- 
haviour Xc^ = ae _A< =/ T , which allows us to extract A c from 
the T dependence of \c- Supplementary material [13] shows 
a plot of ln(xc) vs - The simultaneous strong decrease of 
the low-T Xc an d the opening of A c with increasing U allows 
us to extract the critical value of the interaction U c at which 
the MIT appears. 

In Fig. 2 we see that Xc does not exhibit any sign of diver- 
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Figure 3. (color online) Zero- temperature phase diagram in the U- 
t' plane. Red points with error bars show our estimate of the Mott 
metal-insulator transition (U c ) for various values of frustration t' . 
Above U c is a Mott insulating phase and below U c is a metallic 
phase. Transition at t'/t = corresponds to a square lattice with 
perfect nesting and appears at U c = [21, 22]. For small t'/t, a 
superlinear increase of U c (t') is predicted by the Hartree-Fock ap- 
proximation (shown with black dotted line). 

Phase diagram. In Fig. 3 we show our estimate of the 
critical interaction strength U c (for various t'/t) at which the 
system undergoes a Mott metal-insulator transition. U c de- 
creases with decreasing frustration t'/t and the MIT can there- 
fore be driven either with increasing interaction strength U or 
decreasing frustration (t'/t). At t' = the model is a nearest- 
neighbour square lattice Hubbard model with perfect nesting 
for which U c = [13, 21-23] and going away from perfect 
nesting with increasing t'/t results within a Hartree-Fock ap- 
proximation in a superlinear increase of U c (see Fig. 3). 

Our phase diagram is consistent with previous findings by 
several numerical techniques, e.g., slave bosons and varia- 
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tional Monte Carlo, while it predicts larger values of U c than, 
e.g., path-integral renormalization group, and smaller values 
of U c than, e.g., Brinkman-Rice picture and dynamical mean- 
field theory approaches (see Supp. [13] and Table I therein for 
more details). 

The calculated T-dependence of thermodynamic quanti- 
ties do not show strong signatures of possible different Mott 
insulator spin states (antiferromagnetic order for small if /t 
[24-27], 120 degree Neel order at higher U [25, 28, 29] for 
t' /t ~ 1 and a possible spin liquid at U > U c for t' /t ~ 1 
[26-29]). The spin structure factor or discontinuities in the 
double occupancy [13, 29] would be better indicators [30]. 
Hence, we don't show possible spin ground states in the phase 
diagram. 

Specific heat and entropy. In Fig. 4 we show how the 
T-dependence of the specific heat Cy and entropy per site 
s change with increasing U. In the metallic regime (U < U c ), 
the low-T slope of Cy vs. T and s increases with U. In a 
Fermi liquid picture this slope increase corresponds to the in- 
creased renormalized quasiparticle mass, and we estimate it 
reaches m* /to& ~ 2.5 for U close to U c . m* is the renormal- 
ized quasiparticle mass and mj, is the bare band mass. Com- 
parable enhancements are seen in organic charge transfer salts 
[31]. Simultaneously with increasing slope, the low-T peak 
in Cy moves to lower T, resulting in a decreased coherence 
temperature T CO h- 

We estimate the coherence temperature T m h as the temper- 
ature at which Cy starts to deviate substantially from linear- 
ity in T and obtain T m h < O.li for t' = t in the vicinity 
of the MIT (U ~ U c ). This shows the importance of strong 
correlations, since T co h is much smaller than the estimate of 
T C oh ~ 0.4t for U — 0. Electronic structure calculations 
based on density-functional theory (DFT) give values of t in 
the range 50-70 meV for the k-(BEDT-TTF) 2 X family [32- 
34] and 40-50 meV for the /3'-X[Pd(dmit) 2 ] 2 family [35]. The 
ratio of t' /t varies between about 0.4 and 1.3 depending on 
the counterion X. Taking t ~ 40 meV we estimate T co h < 50 
K, which is in good agreement with experiments [4]. This 
temperature corresponds to the vanishing of the Drude peak 
in the optical conductivity [6], maximum in the thermopower 
vs. temperature [36] or the resistivity becoming comparable 
to the Mott-Ioffe-Regel limit [37]. 

At T > T co h we expect bad metallic behaviour with well 
formed local moments. This is supported by the entropy 
showing already in the metallic regime an increase towards 
the large U result, and furthermore by the large spin suscep- 
tibility close to the result for a Heisenberg model (see Fig. 
5). With further increase of U and entering into the insulating 
regime (U > U c ), both Cy and s are strongly increased at 
low T due to well formed local moments and a large density 
of low lying spin excitations. This is a hallmark of magnetic 
frustration [10]. Closer examination of Cy reveals that the 
low-T peak (at T < O.li) is strongly increased for U = lOi 
and 12t, which might be a signature of a transition from a spin 
liquid state into a Neel order state, in agreement with the find- 
ings in Ref. [29]. High-T properties of Cy and its two peak 
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Figure 4. (color online) Specific heat Cv (top) and entropy s (bot- 
tom) vs. T for several values of interaction U and for t' = t. Top: 
With increasing U, low-T Cv increases. For U <U C this increase is 
due to an increased quasiparticle renormalized mass manifested also 
in an increased slope of Cv at low T and a decreased coherence tem- 
perature T co h. For U > Uc, Cv starts to develop a peak at low T due 
to well-defined low energy spin excitations, while charge excitations 
are gapped with a large charge gap A c , resulting in a high T peak 
(not shown here, see Fig. S3 in Supp. [13]). Bottom: A similar in- 
crease is seen in the entropy, which for U > U c starts to approach the 
value of ln(2). This is characteristic of the development of S — 1/2 
local moments. In particular the strong increase of the entropy at low 
T(~ O.lt) and plateauing below ln2 for T > 0.3t is in contrast to 
what is observed for U = 0: the entropy steadily increases with T 
and only tends towards ln2 at much higher temperatures. Entropy 
is increased at low T even for U < U c , signaling the development 
of local moments already in the metallic regime for U ~ U c , and 
therefore the bad metallic behavior. 

structure at large U (the low-T peak is due to spin excitations 
and the high-T peak due to charge excitations) are shown in 
Fig. S3 of the Supplementary material [13], and for t' = 
were previously discussed in Ref. [12]. 

Spin susceptibility. In Fig. 5 we show, that the spin 
susceptibility \ s is close to the Curie-Weiss (CW) result 
x f w) = 1/(4(T + T cw )) for all calculated values of U, 
in both the metallic and insulating regime close to U c . The 
CW result is obtained by a mean field treatment of a Heisen- 
berg model and T cw = J + J'/2 = At 2 /U + 2t' 2 /U [10]. 
This gives strong support that already for U ~ U c the local 
moment is well formed, resulting in increased \s and a bad 
metal behaviour for U < U c . 
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Figure 5. (color online) Spin susceptibility \ s vs. T for several val- 
ues of U and t'/t = 1. \ c increases with increasing U towards 
the Curie law Xs(U = oo) = 1/(4T) or noninteracting spin 1/2 
result. For intermediate U, Xs is close to the Curie- Weiss (CW) 

law xf W) = 1/ 4 ( T + Tew) with Tew = J + T/2 (shown 
for U — 12t and 5i) and therefore the behaviour is close to the mean 
field result for the Heisenberg model. Deviations of Xs from the CW 
result appear at high T > 1.5t due to increased double occupancy 
and decreased local moment, while at low T < 0.5t deviations occur 
for U > U c due to increased longer range correlations. The strong 
decrease of Xs at T < O.lt and U > U c could also be due to strong 
increase of spin correlations or the opening of a gap to triplet ex- 
citations. Xs is close to the CW result even in the metallic regime 
for U < Uc and it does not show the weak T dependence charac- 
teristic of a Pauli susceptibility found in Fermi liquids, consistent 
with our estimate T co h < O.lt. These results suggest well formed 
local moments with short-range antiferromagnetic correlations and 
bad-metallic behavior. Unlike the charge susceptibility Xc the spin 
susceptibility Xs does not show any signature of the metal-insulator 
transition. 
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Figure 6. (color online) Effect of frustration (t'/t) on thermodynamic 
properties in the Mott insulating phase. Here we show the tempera- 
ture dependence of the specific heat CV (top), entropy s (bottom left) 
and spin susceptibility Xs (bottom right) for U = lOi > U c and sev- 
eral values of t'/t. With increasing frustration (t'/t) the low energy 
density of states is strongly increased as seen in the strong increase 
of C v and s at T < OAt. These states are low energy spin exci- 
tations which can be seen in increasing of s towards In 2 (free spin 
result) and increased low-T Xs ■ Dramatic shift of the peak in Cv 
and Xs towards lower T with increasing frustration can be attributed 
to decreased spin correlations. This and the increased density of low 
energy spin excitations is a hallmark of magnetic frustration [10]. In- 
creasing t'/t beyond 1 starts to reduce the frustration and the peak in 
Xs again moves to higher T (shown for t'/t — 1.5). 



Agreement with the CW behaviour supports the picture of 
a bad metal with short range antiferromagnetic correlations, 
while longer range correlations may suppress Xs below CW 
at quite low T. Correlations start to develop below 0.5i and 
become stronger for T < O.lt. Alternatively, the suppression 
of Xs for T < O.lt might occur due to a spin gap to triplet ex- 
citations. Such suppression of longer range spin correlations 
is due to magnetic frustration [10]. 

Even in the metallic phase we observe a suppressed Xs at 
low T (see Fig. S7 in Supp. [13]). A similar suppression has 
been observed in the T dependence of the NMR Knight shift 
in some organic charge transfer salts [38]. It has been argued 
that this suppression together with a similar suppression in 
the NMR relaxation time 1/T X T [38-40] is a signature of a 
pseudogap. We further discuss this issue in the Supp. [13]. 

In Fig. 6 we show how frustration affects the low-T proper- 
ties of the insulating phase (U > U c ). Frustration strongly in- 
creases the density of low-lying spin excitations, which results 
in a dramatic increase of Cy and s at low T. Furthermore, Xs 
is also increased at low T due to suppression of longer range 
correlations and its peak moves to lower T. This is in agree- 



ment with series expansion results [10] for the corresponding 
Heisenberg models. Such dramatic increase of the low-T Cy, 
s and Xs is a hallmark of magnetic frustration. In the bottom- 
right of Fig. 6 we also show, how with increasing t' beyond t 
and subsequently decreased geometric frustration, the peak in 
Xs moves to higher T again. 

In conclusion, we have considered the thermodynamic 
properties of a Mott MIT, which can be driven either by inter- 
actions (U/t) or geometric frustration (t'/t). We have shown 
that the metallic phase near the MIT is characterized by a 
small charge susceptibility, large quasiparticle renormaliza- 
tion, a reduced coherence temperature T CO h ~ O.lt, large en- 
tropy, and large spin susceptibility. This is in agreement with 
experiments on organic charge transfer salts [4]. We have ar- 
gued that the large entropy is due to a large local fluctuat- 
ing magnetic moment, which leads to bad metallic behaviour 
above T co h. Furthermore, we have shown how frustration in- 
creases the density of low energy spin excitations and reduces 
the range of antiferromagnetic spin correlations in the insulat- 
ing phase. 
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FINITE TEMPERATURE LANCZOS METHOD 

In this section we present in more detail the finite tempera- 
ture Lanczos method (FTLM) [1,2] and values of some of the 
parameters used. In particular, we focus on the accuracy of 
the method, for which several parameters need to be consid- 
ered. These are the number of Lanczos states M, the number 
of samples over starting random vectors R, and the system 
size N. 

The number of Lanczos states M used in our calculation 
varied from 50 to 100, which is orders of magnitude smaller 
than the number of basis states of a 16 site cluster (in which 
the Lanczos states are written). Such Ms are sufficient for ob- 
taining convergence of the result. The ground state converges 
within such M as well as finite T properties. This can be 
traced back to the fact that moments up to the order of M are 
exact for a state of interest (see eq. 3.18 and corresponding 
text in Ref. [1]). Thermodynamics is only weakly dependent 
on M as can be seen in top Fig. 3 in Ref. [1], where already 
M — 5 and 20 gave quite accurate results. When dealing with 
dynamics and spectral properties, M limits the frequency res- 
olution and usually larger values of M are employed (Fig. 4 
and 5 in Ref. [1]). 

In contrast to the zero T Lanczos method, one employs 
within FTLM averaging over random vectors in order to cal- 
culate finite T properties. This is most nicely described and 
justified in section 3.5 in Ref. [1], It is shown that the relative 
statistical error is of the order of 



0(- 
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where R is the number of random vectors used and 

Z = TrcM-P(H -E )). 



(SI) 



(S2) 



Z is the thermodynamic sum normalized by exp(— (3E ), 
where Eq is the ground state energy. Therefore the error is 
very small for large Z, which can appear either at high T or 
for larger systems. Using larger systems with larger Z reduces 
the error as well as finite T, since Z can strongly (e.g. expo- 
nentially) increase with increasing T. Errors usually become 
larger at low T, where one would need to employ a large R 
to reduce the error. In our case averaging over random vec- 
tors was in combination with averaging over twisted bound- 
ary conditions with 32 or more samples. Also, as shown in 
Fig. 6, frustration increases the density of states and entropy 
at low T and helps improve accuracy in this respect. Some 
dependence on R can be found in the bottom panel of Fig. 3 



in Ref. [1], where a system with smaller Hilbert space (and 
therefore smaller Z) was used. Therefore the most computa- 
tionally demanding regime within FTLM is low T, which can 
be made less demanding by improvements suggested with the 
low temperature Lanczos method (LTLM) [3]. 

More challenging are the finite size effects, which are 
largest at T = (potentially long correlation lengths), but 
become smaller with increasing T, since correlation lengths 
decrease with increasing T. In Ref. [1] (section 3.7) it is ar- 
gued, that finite size effects are small for T above Tf s , at which 
Z reaches a certain value (<~ 30). At such an elevated T, sys- 
tem size (7Y) dependence of the results becomes small and 
one essentially obtains a result close to the result for an infi- 
nite system. We show some system size dependence in Figs. 
SI and S2. Frustration (present in our model) reduces the spin 
correlation length and also reduces the relevant energy scales 
(for example see Fig. 6), which makes Tf S smaller. Further 
reduction of the finite size effects can be obtained with aver- 
aging over twisted boundary conditions [4, 5] (Fig. SI). 

To shortly summarize the above discussion one can use the 
approximation for (S2) 



Z(T) ~ exp[Ns(T)], 



(S3) 



where N is the number of sites and s(T) is the entropy per 
site. From this we see, that at fixed (low) temperature s(T) 
is significantly increased by the interactions (Fig. 4) and the 
frustration (Fig. 6) and so this (i) reduces finite size effects, 
(ii) extends to lower temperatures the regime of reliability of 
the numerical method, and (iii) reduces the statistical noise. 
Thus, the numerical method is most reliable in the parameter 
regime of greatest physical interest: strong interactions and 
large frustration. 



TESTING THE NUMERICAL METHOD 

In Fig. SI we show how averaging over twisted boundary 
conditions improves results and makes the finite size effects 
smaller. In Fig. S2 we show the system size dependence of 
the specific heat Cy. 



FURTHER DISCUSSION OF X c 

At zero temperature and for U = 0, \c = 2/Vo(/i), where 
No(n) is the non-interacting electron density of states. In Fig. 
1, Xs for the non-interacting electron case (U — 0) shows a 
peak at T ~ 0.4t due to the van Hove singularity. 
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Figure SI. (color online) Averaging over twisted boundary condi- 
tions reduces finite size effects. The temperature dependence of the 
charge susceptibility is shown for different numbers of samples over 
boundary conditions Ng for U = 6t and lOt. For U — 6t and 
Ng = 1 and system size N = 16, \ s shows a strong increase for 
T < QAt, which is a finite size effect. For example, N = 14 results 
for large Ng shows no such increase and actually shows suppression 
for T < 0.2t. With increasing N$, \c for N = 16 no longer shows 
strong increase at low T and becomes only weakly T dependent as 
expected for a metal and is therefore closer to the result in the ther- 
modynamic limit. Changes of \c with larger Ng are small in the 
insulating phase (U = lOt), which can be traced back to the van- 
ishing of the Drude weight, since the Drude weight can be calculated 
with the derivative of the energy with respect to the twisted boundary 
phase 9 [6-8]. 



Xc decreases with increasing U, which is most apparent at 
high T (~ 2t) due to spreading of the density of states over the 
larger energy interval of W + U or i.e. over a larger chemical 
potential interval. Here W is a bare band width. The chemical 
potential interval \p(n = 0),/x(n = 2)] in which density of 
electrons changes from n = to n = 2 can be determined 
exactly. /x(n = 0) = E (l) - E o (0) = ef in , where E (N e ) 
is the ground state energy of a system with iV e electrons and 
£ mm j s tne ener gy a ( the minimum of the bare band. Similarly, 
H(n = 2) = E (2N) - E (2N - 1), where N is the system 



size, E {2N) = NU and E (2N - 1) = NU - e£ 



U. Eq(2N — 1) can be easily calculated with the particle- 
hole transformation a. a <H> c- a , 
with t 



T , resulting in the Hamiltonian 
tij and an extra term U(N — N e ) [9] and for 
Eq(2N — 1) only one particle state needs to be considered. 
/Lt(n = 2) = t^} ax + U . Therefore the interval of the chemical 
potential in which n rises from to 2 is increased by U to 



11" 



U, resulting in on average decreased \ c — At low 
T and metallic regime, Xc can be related to the quasiparticle 
weight z, which is discussed in the next section, while at high 
T, X c(T » W) = 1/(2T). 

In Fig. S3 we show ln(x c ) vs. 1/T, which makes the 
opening of the charge gap A c in the insulating phase clearly 
seen. Whether linearity of A c (f7) extends all the way to U c 
or A c = cannot be concluded from our FTLM results, since 
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Figure S2. (color online) System size dependence of the specific heat 
Cv and entropy s. For systems with iV ~ 16, Cv shows only small 
quantitative changes with system size at intermediate T(~ 0.3t), 
while at high- and low-T shows only weak dependence on the system 
size. 



small gaps affecting low temperatures cannot be reliably ex- 
tracted due to the finite size effects. Our result of the opening 
of A c linearly with increasing U is also consistent with T = 
exact diagonalization results in Ref. 10. Although their values 
of A c are larger than ours by approximately 0.3t, they show 
a trend towards our values with increasing system size. An 
approximately linear in U opening of the gap was also found 
for a Hubbard model on a square lattice with n next nearest 
neighbor hopping t' Jt = 0.2 using the grand-canonical path- 
integral renormalization group [11]. 



QUASIPARTICLE RENORMALIZATION AND CHARGE 
SUSCEPTIBILITY 



Here we consider the effect of a quasiparticle renormal- 
ization on the charge susceptibility within a simple model 
and show that for a simple Fermi liquid, Xc is proportional 
to a quasiparticle weight z. This is in contrast to what is 
naively expected from the increased low energy density of 
states due to increased renormalization or quasiparticle ener- 
gies (e fe -> ze k ), e.g. Xc oc 1/z. 

The density of electrons n in a system with N sites can be 
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Figure S3, (color online) Arrhenius plot of ln(Xc) vs. 1/T for sev- 
eral interaction strengths U and t' — t. In the insulating phase the 
activated behaviour is linearized in this plot, ln(xc) = ln(a) — A c ^ , 
with the slope of the curves given by — A c . This makes the opening 
of the charge gap A c nicely seen as increased negative slope of the 
curves. 



calculated as 



2 

N 



k 



(S4) 



The factor 2 is due to spin, while n k can be obtained from 
the spectral function or imaginary part of the retarded Green's 
function, 



rik = / du>f(oj)lm — 

7T ./ uj - 



1 



[i - e k - S fc (cj) ' 



(S5) 



f(oj) is a Fermi-Dirac distribution function, f(uj) — 1/ (e^ + 
1), 6^ is a bare -band dispersion, fi is a chemical potential, 
and E/j (lu) is a self-energy. Using the definition of the charge 
susceptibility Xc 



9jU 



one can write 



x - dn k 
~N ~dj~L 



(S6) 



Furthermore, can be expressed in terms of the real S' fc (w) 



and imaginary S' fe '(w) parts of the self-energy, 



On 



/(w)E£(w)(w + /i-e fc -E / fc (w)) 



(S7) 

Here we have assumed that the self-energy has no depen- 
dence on ii. Since we are interested mostly on effects of 
renormalization on Xc within the quasiparticle picture we also 
use the following approximations. We approximate T^ k (uS) = 
9 W £' (lu^^qlu, which leads to 



u + fj, — ek — S'(w) ~ — + fi - e k , 
z 



(S8) 



where z = 1/(1 — 9 w E'(u;)| u= o) and assume S^w) ~ S" 
or, with other words, we neglect the to and k dependence of 



E' fe '(w) close to the quasiparticle peak. With this approxima- 
tions we can perform the integral in Eq. (S7) in the T — > 
limit. 



— (i — > 0) = Z 7T. 



(S9) 



In a quasiparticle picture we assume | S" | -C |ej. — which is 
satisfied in the most relevant regime close to the Fermi surface 
and leads to the approximation 



^(T^0)~z<5(e fe - M ). 



(S10) 



The effect of quasiparticle renormalization on \c now be- 
comes clear, 



Xc(r^0) = |^z<5(e fe - A1 )^zx! 

k 



(Sll) 



where Xc = 2JVo(m) ' s tri e bare charge susceptibility. Charge 
susceptibility is due to quasiparticle renormalization reduced 
from it non-interacting value by a quasiparticle weight z. 

This simple model shows that due to interactions and quasi- 
particle renormalization Xc is reduced by a factor z, which is 
what we qualitatively observe in our numerical results (Figs. 
1 and 2 in the main text) where we assume that at low temper- 
atures CV ~ T/z. This is in agreement with DMFT results 
[12], but conflicts with a claim in Ref. [13], where they sug- 
gest that both Xc and the specific heat coefficient 7 should be 
proportional to l/z. This simple model and our results are 
consistent with Xc oc z and 7 oc l/z. 

The divergence of Xc with approaching a Mott insulator by 
reduced filling [13-15] can not be captured with this simple 
model, and might be due to strong dependence of the self- 
energy on fj, and/or breakdown of a quasiparticle picture. 



FURTHER DISCUSSION OF THE PHASE DIAGRAM 

Deviations from perfect nesting and t' /t = results in a 
strong exponential or superlinear increase of U c . In Ref. [16], 
this was observed by shifting the chemical potential and the 
dependence of the critical chemical potential fi c on U is given 
by /i c ~ VW exp(— 2iT\ft/U). Furthermore, strong increase 
of U c with increasing next nearest neighbor hopping in the 
square lattice was also observed within a Hartree-Fock ap- 
proach [17] and a similar result is obtained for small t'/t on 
the anisotropic triangular lattice as is shown in Fig. 3 with 
black dotted line. 

Describing the behaviour of U c (t') at small t' seems more 
challenging, since Hartree-Fock approximation gives a super- 
linear increase of U c with t' [17], exact diagonalization [18] 
and VMC [15] suggest a linear increase, while cellular DMFT 
[19] result seems to be more consistent with quadratic-t' de- 
pendence. 

In Table I we compare critical values of U c for MIT as ob- 
tained by different methods. DMFT denotes dynamical mean- 
field theory. 
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Table I. Comparison of the critical interaction strength U c for a MIT 
ait' — t as estimated by different methods. 



Method Uc/t Reference 



Slave rotors 




5.14 


[20] 


Path-integral renormalization group 




5.2 


[21] 


Hartree-Fock 




5.27 


[22] 


Strong coupling expansion 




6.7 


[23] 


Variational cluster approximation 




6.7 


[24] 


Exact diagonalization for T = (N = 


16) 


7 


[18] 


Slave boson with magnetic order 




1 .2.5, /.Do 


|25, lb] 


FTLM 




7.5±0.5 


this work 


Variational Monte Carlo 




7.65±0.05 


[27] 


Cellular DMFT 




8.5, 10.5 


[19, 28] 


Cluster DMFT 




9.2-9.6 


[29, 30] 


DMFT 




12-15 


[31,32] 


Exact diagonalization for T = (N = 


12) 


12 


[26] 


Resonating-valence-bond theory 




12.4 


[33] 


Brinkman-Rice 




15.8 


[26] 
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Figure S4. (color online) Specific heat CV vs. T for several values 
of interaction U and for t' = t. With increasing U, Cv develops 
two peaks, one at low T, which is at large U > U c due to spin 
excitations, and one at high T corresponding to charge excitations 
across the charge gap to the upper Hubbard band. 

For U ~ U c , Cv already starts to develop a two peak struc- 
ture (see Fig. S4), which is most pronounced for U deep in 
the insulating regime (U ~ 12t > U c ). The low T peak 
corresponds to spin excitations as was already observed for a 
Hubbard model on a square lattice [5], where the low-T peak 
appeared at T ~ J and for large U can be captured with the 



Heisenberg model. In our case the peak is moved to lower 
T (< J) due to frustration. The high-T peak corresponds 
to charge excitations across the charge gap A c into the upper 
Hubbard band. 

DOUBLE OCCUPANCY AND LOCAL MAGNETIC 
MOMENT 

The double occupancy D can be calculated with the use of 
the free energy F defined by 

e-P F =Tre-P H . (S12) 

Usually one evaluates D by taking the derivative of F with 
respect to U at fixed chemical potential //. On the other hand, 
we are dealing with a fixed number of particles or fixed filling 
and therefore in our case fi changes with U. Taking this into 
account, one gets 

1 OF d ll 

Therefore, D is calculated from the derivative of F with re- 
spect to U at fixed number of electrons N e , and one needs 
to add a term due to the change of chemical potential with U 
at fixed N e . N is a number of sites in the system. Our cal- 
culation of D serves only as a rough estimate, since we take 
numerical derivatives of F and /i for quite large AU ~ t. This 
does not allow for a precise determination of D, and smooths 
out any sharp features of D as a function of U. 

In Fig. S6 we show the calculated U dependence of a lo- 
cal moment (si), which shows a smooth behaviour without 
any substantial change near U c . This supports the picture of a 
large local moment, with values close to the strong coupling 
(Heisenberg) limit even in the metallic phase. Our results do 
not show the behavior of a local moment predicted with the 
Hartree-Fock or Slater approximation, where the local mo- 
ment is a constant with the non-interacting value for U < U c , 
and increases slowly with increasing U for U > U c . In 
this approximation the MIT is driven by antiferromagnetism, 
while our results are consistent with the MIT driven by Mott 
physics. Our results are also in contrast with the Brinkman- 
Rice picture [34], which predicts that the local moment in- 
creases linearly with increasing U for U < U c ~ 15.8t, 
fully develops (4(s^) = 1) at U — U c and stays constant 
for U >U C . 

The agreement with the strong coupling result seems sur- 
prisingly good for U close to U c , which suggests that the 
Heisenberg model gives a good approximation also in the 
regime U > U c and that higher order terms do not play a 
crucial role. This appears in contrast with results in Ref. 23, 
where they observed that higher order terms are actually re- 
sponsible for the transition between Neel ordered and spin 
liquid phase. However, our results are for finite T where the 
small differences in the ground state energies are not that im- 
portant, and also the change of D at the transition was ob- 
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Figure S5. (color online) Double occupancy D vs. interaction 
strength U for t' = t and T — O.lt. D decreases with increas- 
ing U as expected. For small U < U c , D shows linear decrease with 
increasing U, which is predicted with the Brinkman-Rice (BR) pic- 
ture [34] (D (BR) = (1 - U/U^ BR) )/4:). However, the BR picture 
predicts too large a Ui BRS> ~ 15.8i [26, 34]. Close to U c , D is more 
strongly suppressed and starts to exhibit strong coupling (SC) [9, 35] 
or Heisenberg behaviour in which double occupancy is given by [35] 
D (SC) = (2t 2 /U 2 ) Eit 1 / 4 - < s i ■ Si+t})- The sum over 6 goes 
over all 6 nearest neighbours and the strong coupling result is shown 
for (Si ■ Sj+j) ~ —0.182 [36]. This estimate of spin correlation is 
evaluated within Heisenberg model and is valid for U > U c , where 
it shows only a weak dependence on t/U [9] . The agreement of 
the calculated D and the SC result is surprisingly good in the regime 
shown in the figure (U > U c ). The small value of D ~ 0.1 close to 
U c corresponds to only every tenth site being doubly occupied, which 
results in a large local moment and strong spin response manifested 
in large \a ■ We note that in Ref. 23 a small discontinuity (~ 0.01) in 
D was observed at U ~ Wt and attributed to a first-order transition 
from a spin liquid (6t < U < Wt) to a Neel antiferromagnet with 
120 degree spiral order (U > Wt). Our results do not have sufficient 
resolution to detect such a transition. 

served to be only a few percent [23]. Furthermore, we esti- 
mate D by numerically differentiating the free energy over U 
with AZ7 ~ t, which further smooths the [/-dependence of D. 

PSEUDOGAP 

An important question is whether a pseudogap is present 
in the metallic phase close to the Mott insulator [38]. Sig- 
natures of such a pseudogap are seen in NMR experiments 
on K-(BEDT-TTF) 2 Cu[N(CN) 2 ]Br [37, 39]. The Knight shift, 
which is proportional to the spin susceptibility Xs> decreases 
by about 40% as the temperature is lowered from about 50 K 
to 10 K. The NMR relaxation rate 1/T±T increases with de- 
creasing temperature, with a maximum around 50 K, and then 
decreases by about a factor of two as the temperature is low- 
ered to 10 K. These temperature dependences are qualitatively 
similar to what is observed in underdoped cuprate supercon- 
ductors for which ARPES provides independent evidence of a 



Figure S6. (color online) Local moment 4(s^) = 1 — 2D vs. interac- 
tion strength U for t' /t — 1 and T — O.lt. l{s 2 z ) is increasing with 
increasing U as expected and has values close to the ones expected 
from the strong coupling (SC) limit [35] for U > U c . Large values 
of the local moment persist also in the metallic regime for U < U c 
and we do not observe a strong decrease of (s 2 ) with decreasing U 
at the MIT. Therefore the metallic phase is characterized also with 
a large local moment and strong spin response, e.g., with large \s- 
The Brinkman-Rice result [34] is also shown, together with limiting 
values for U = oo and U = 0. 

pseudogap. 

Figure S7 shows that for t' = 0.8t and U = 6t, Xs de- 
creases by about 50% as the temperature decreases from about 
O.lt to 0.03i. These parameters correspond to the metallic 
phase, as indicated by the non-zero charge compressibility. 
The calculated temperature dependence appears to be consis- 
tent with the experiment [t ~ 500 K]. However, caution is in 
order, because of the possible role of finite size effects at such 
low temperatures. This can be seen by examining the tem- 
perature dependence of Xs for U — lOt which is in the Mott 
insulating phase. It has a maximum around T = 0.15t and 
decreases smoothly to zero around 0.03i. Similar behaviour 
is found for U — lOt and smaller values of t' /t. These re- 
sults can be compared to known results for the corresponding 
Heisenberg model. In particular for J' < 0.5J [t' < 0.7t] 
the model should have long-range Neel order at zero temper- 
ature. Quantum Monte Carlo simulations on the square lattice 
model show the temperature dependent spin susceptibility has 
a maximum around T ~ J and then decreases by about 50% 
to a non-zero value at T — 0; hence, there is no spin gap 
[40]. This discrepancy shows that the apparent gap observed 
in FTLM is a finite size effect. On the other hand, the sup- 
pression of Xs, m the temperature range 0.05i < T < O.lt 
may be a real effect. But, it is not clear at a moment what is 
the physics behind this reduction in the spin susceptibility. It 
could be due to a suppression of the density of states such as 
associated with a pseudogap. Or like in the Heisenberg model 
the reduction could be due to the development of longer-range 
antiferromagnetic correlations in the bad metallic phase. 
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Figure S7. (color online) Sign of a pseudogap in the metallic phase 
close to the Mott MIT. The temperature dependence of the spin (Xs) 
and charge susceptibilities (Xc) are shown for t'/t = 0.8. The sup- 
pression of Xs for U — 6t and 0.03t < T < O.lt (top panel) could 
be due to the emergence of a pseudogap. In this regime the system 
is metallic, as can be seen in the T independent Xc (bottom panel). 
This is in agreement with the measured Knight shift K s [37], How- 
ever, caution is in order, since the suppression of Xs for U = lOt at 
T < 0.05i may be over-estimated due to finite size effects (see text). 
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